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ABSTRACT 


This report summarizes the analysis and results developed in a fresh 
approach to calculate flow induced vibration response of a flexible flow 
passage. The vibration results are further examined in the frequency domain 
to obtain dominant frequency information. A cumulative damage analysis due to 
cyclic strains is performed to obtain number of cycles to failure for metallic 
bellows of particular specifications under a variety of operational 
conditions. Sample plots of time and frequency domain responses are included 
in Appendix I. Complete listing of a computer program is provided in this 
report as Appendix II. The program successively executes each of the analyses 
needed to calculate the vibration response, the frequency response, the cyclic 
strains and the number of cycles to failure. The program prompts the user for 
necessary input information. Sample data from the program is provided in 
Appendix III. The fatigue life results obtained by the computer model lie 
within an acceptable range of previously measured available data. 



I. INTRODUCTION 


Project Background and Overview 

This report describes all the work performed by Georgia Institute of 
Technology, College of Engineering, School of Mechanical Engineering under 
Grant number NAG10-0017, "Fatigue Behavior of Flexhoses & Bellows Due to Flow 
Induced Vibrations". This study was performed for the Kennedy Space Center of 
the National Aeronautics and Space Administration. 

The general objective of this work was to investigate the physical 
phenomena associated with flow-induced vibrations of a flexible line in order 
to develop a methodology to examine its fatigue life. In response to NASA 
needs, much emphasis was placed on developing a computational procedure for 
calculating the number of cycles to failure for a specific bellows based on 
the assumed vibration model. 

Flexible expansion joints are currently being used in a variety of 
applications. These include the nuclear reactor cooling systems, supply and 
return loading lines for submarines from the mother ship on the base, liquid 
cryogenic fuel lines for the space shuttle external tanks and so forth. There 
is, in fact, a manufacturer's association for expansion joints which publishes 
its recommended practices for their use. 

Serious difficulties are encountered when practical solutions are not 
obtained to control the so-called "self-generated" sound or vibrations in flow 
systems which use such corrugated lines. For the special case of the 
internally corrugated flexible hoses, segements of which are utilized in the 
feed lines of Space Shuttle external tanks, a qualitative mechanism was 
identified in terms of the fluid flow related parameters governing the 
physical phenomena [1], Unshrouded hose vibration was attributed to an 
approximate matching between the longitudinal structural frequency and the 
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frequency of vortex shedding at a critical flow velocity over the first few 
corrugations in the hose. This led to the computation of the modal 
frequencies of vibration by representing the hose with a series of lumped 
spring-mass elements and relating the amplitude of vibrations to the damping 
characteri sties of the system. It was found necessary to supplement the 
postulated model by a large number of empirical correlations obtained from a 
series of experiments to relate the system geometry, the mean flow rate and 
the fluid properties with the modal frequencies and amplitudes of vibration, 
the damping characteristics, the maximum alternating stress and the expected 
fatigue life of the hose material [2], 

Experimental observations under simulated field conditions have proved 
the predicted failure results of the two cited works to be very 
conservative. Furthermore, an overwhelming reliance on empirically obtained 
functions to predict the fagitue life of the hose in these studies has made it 
cumbersome to directly utilize the suggested failure prediction procedure. 
Flex line flow-induced vibrations were also examined in another study [3] with 
an eye to developing a methodology to assess its fatigue life. The cavity 
resonance model utilized in this work tends to be extra sensitive to minor 
parameter variations in predicting the time to failure by fatigue stresses. 
It appears that a definite need exists in developing a vibrational model for 
the flexible line that incorporates the interactive dynamics between fluid 
flow past cavities and the flow-induced vibration of the convolutes of the 
bellows configuration. 

In the context of flow past cavities, there are three types of self 
generated oscillations. These are the fluid-dynamic, the fluid-resonant and 
the fluid elastic types of oscillations [4,5]. The fluid dynamic oscillations 
arise from some type of flow instability inherent to the flow configuration. 


- 3 - 



The fluid-resonant oscillations occur due to standing wave resonance effects 
which are particularly important for compressible fluid flow. The fluid- 
elastic oscillations involve a coupling of the solid boundary vibrations with 
the fluid flow perturbations causing those. It is this coupling which must be 
quantitatively included in the vibration model which predicts the convolute 
deformations. Subsequent cyclic strain analysis and the fatigue damage due to 
it are also examined in this work. 

Statement of Research Problem 

This research concerns the development of a vibration model as well as 
fatigue damage analysis to estimate life cycles of a bellows or a fl exhose 
with a single- or multiply wall construction. The vibration model emphasises 
the coupled interactive nature of the fluid flow and structural vibrations. 
The vibration response is further examined to assess the dominant frequencies 
causing fatigue damage. Finally, the research adopts existing theories of 
cumulative fatigue damage with the results of vibration analysis to assess the 
number of cycles to failure of the bellows. 

Flow Past Cavities 

One of the first modifications to improve currently existing models of 
flow excitation in a bellows geometry concerns the force acting on each 
individual convolute. A preponderence of previously published information 
concerns flow past rectangular cavities. Of particular interest is the 
pressure distribution along the walls of the cavity. 

Sketches of flow patterns in a rectangular cavity based on photographs 
[6] are shown in Figures 1(a), (b) and (c). 

A plot of the measured pressure distribution [6] along the cavity wall 
for three different values of aspect ratio, h/b (0.5, 2.0, 1.5) reveals a 
trend shown in Figure 2. Similar plot [7] is shown in Figure 3. On the other 
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Figure 1. Sketches of Flow Past Square Cavities Based on 
Photography [6]. 
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Figure 2. Typical Measured Pressure Distributions Along Cavity 
Wall Showing Influence of Height-to-width Ratio [6']. 
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Figure 3. Pressure Distribution Along Square Cavity Walls [7] . 
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Figure 4. Pressure Distribution From Analytical Model [6J. 
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hand, an analytical model and its results of pressure variation along 
horizontal axis [6] are shown in Figure 4. It is seen that there is a 
pressure maximum at an area of flow reattachment on the downstream wall of the 
cavity. 

On the assumption that these same trends hold for flow past a convolute- 
shaped cavity, the expected flow patterns are as shown in Figure 5. The 
maximum wall presure is shown acting near the top of the downstream wall as 
demonstrated in the earlier results on cavity flow. The drop in pressure just 
previous to the sharp pressure rise shown in the plots might be explained by 
the vortices in the recirculating flow on the downstream wall. The vibrations 
of the bellows cause these reci rculating vortices to move along the wall and 
be shed from the convolute tip periodically. As a vortex crosses the region 
of maximum pressure it causes a drop in this pressure. After the vortex has 
been shed the maximum pressure is reestablished, thus producing a fluctuating 
force. This fluctuating pressure force, along with an acoustical resonance 
(radial) that may exist, could produce relatively high pressure levels 
resulting in a bending of the convolute. 

The maximum pressure force pictured may be broken down into transverse 
and axial components. Previous research has described the force acting on the 
convolute as a longitudinal force. Actually, this longitudinal force is seen 
in the present research as only the horizontal component of the maximum 
pressure force. 


II. DEVELOPMENT OF THE MATHEMATICAL MODEL 

Assessment of fatigue failure of metallic bellows due to flow induced 
vibrations requires 

a. the development of a coupled fl u id- structure interactive vibration 
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Figure 5. Impact Pressure Due to Flow Past a Convoluted Cavity. 
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model to predict individual convolute deformations, 

b. a frequency domain analysis of the vibrational response to identify 
dominant frequencies, causing the damage 

and 

c. a failure analysis based on available cummulative damage theories to 
calculate the number of cycles to failure. 

Flow Induced Vibration Model 

The interactive dynamics of bellows vibrations due to fluid flow over a 
cascade of convoluted cavities, which in turn modifies the fluid flow and 
forces that excite the vibrations, is a coupled nonlinear phenomena. It is 
governed by fluid flow conditions, the bellows geometry and the material 
response behavior. 

The effective excitation force is internally generated and an analysis 
must be performed to express it in terms of identifiable physical 
quantities. The vibration of an individual convolute within a multi convoluted 
bellows in response to the excitation may be examined via lumped springs and 
masses in terms of physical quantities. 

The Excitation Force 

A comprehensive analysis of flow induced vibration of a bellows must 
consider the motion of an individual convolute in two coordinate directions. 
These directions correspond to the longitudinal and the transverse motion of 
the convolute as illustrated in Figure 6(a). 

An expression for the excitation force acting on the convolute must 
account for the variation of the force due to convolute motion. When the 
bellows are absolutely motionless the flow geometry is that illustrated in 
Figure 6(b). This flow may be represented by the velocity triangle in Figure 
7(a). 

In this figure, a Q , estimated from data on flow past cylinders [8],is 
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Figure 6a. Coordinate System to. Figure 6b. Configuration for Unde-- 
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Figure 6c. Convolute Deformation Along x. 



Figure 6d. Convolute Deformation Along y. 
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the angle between the velocity of the flow approaching the convolute tip, U 0 , 
and the velocity of the core flow, U c . The flow configuration illustrated in 
Figure 6(b) produces force coefficients on the convolute tip given by 

u 2 

C 0 ,x " p- < C D cos »o + C L s, % ) (la) 

and C 9 

V 

C o,y = TTT ^ C D sina o + C L COSa o^ » 
u c 

where C D and Cj_ are drag and lift coefficients, respectively, of convolute tip 
in the core flow; C D and C L are estimated from flow past cylinders [8]. 

When the convolute tip is displaced longitudinally, the configuration 
illustrated in Figure 6(c) results. The angle made by the flow approaching 
the convolute tip with the core velocity, (* a Q + Aa) is shown on the 
velocity triangle in Figure 7(b) representing this situation. In this 
figure, U -is the velocity of the flow approaching the i th convolute tip 
relative to the convolute longitudinal velocity. The force coefficients for 
this situation become 


c *,i 


_ U r,x,i 


U 


( C D cosa^ + C L si net..) 


and 


Cy i = r> ^ - (C R sin a i + C L cos a.) 


(2a) 

(2b) 


These expressions show that longitudinal motion of the convolute produces a 
variation in both the longitudinal and transverse excitation forces. 
Therefore, longitudinal motion cannot exist independent of transverse motion . 
When the convolute tip is displaced in the transverse direction, the 
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Figure 7a. Velocity Triangle: Un- 

deformed Bellows. 




Figure 7c. Velocity Triangle: Deformation Along y. 



Figure 8. Velocity Triangle Combined x & y Deformation. 
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configuration shown in Figure 6(d) represents the situation. The new 
angle, 8^, made by the approaching flow with the core velocity (8^ = a Q + A8) 
is shown on the velocity triangle in Figure 7(c). In this figure U r ,y,i is 
the velocity of the flow approaching the i t * 1 convolute tip relative to the 
convolute transverse velocity. The force coefficients for this variation are 


and 


If 

C x,1 ' “Hr (C L SlnS 1 +C D C0S «i> 


u 


u 


C y>i * r - ’-y (C L C0S8 i + CpSine^ 


U 


(3a) 

(3b) 


Similar to the case for longitudinal motion, transverse motion causes 
variation in both the longitudinal and transverse forces acting on the 
convolute tip. Therefore, transverse motion cannot exist independent of 
longitudinal motion . 

If the transverse and longitudinal displacements are allowed to occur 
simultaneously, the velocity triangle becomes that shown in Figure 8. From 
this figure the angle <)>.. and the relative velocity U r j may, respectively, be 
determined as 


U tana 0 + (y, - y,^) 

, = tan L J 

U c - <*,_! - x t ) 


(4) 


and 


11^. = [U c tanu Q + (y^ ~ ^i-1^ + ^c ” ^ x i-l - 


(5) 


The force coefficients are given by 



2 ' si n<j>^ + Cp co s<^^ 1 
C 


(6a) 
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and 


C . 

y,i 


u 2 . 

r,i 


2 ~ [ C L cos$. + C D sinfy] . 


The excitation forces may, therefore, be expressed as 


and 


F x,i ■ 7 A C x,1 


F y,1 ‘ W*y.1 ’ 


(6b) 


(7a) 

(7b) 


where pis the fluid density and A is the convolute tip surface area. 

Upon substituion of the expression for <j>^ into the expressions for the 
force coefficients one obtains 


_ r,i 


'x,i 


{ C L sin [tan' ( 


1 , u c tan “o + <*1 - Vl> 


U C ' ^*1-1 " 

-1 , u c <*i “ >1-1* 


)] 


+ C n cos[tan~ ( 

D U. - (x 


1-1 -* 1 > 


) 3 ) 


and 


U 2 . 
_ r,i 


, U r tana + (y - y .) 
C i = ~ L T~ { C L cos Ctan ^ ; : 

.1, « c tan “o + ' *i-i ! 


+ Cp sin[tarf ( 


U C ' <*1-1 ' *i } 


)]} . 


(8a) 


(8b) 


Upon introducing the trigonometric identities 


and 


sin(tan' 1 e) 


e 

l + e 


2 
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cosltarf^e) - 1 


1 + e* 


the force coefficients, with some rearranging, become 


vt . u.tanct + (y, - y, ,) 

■ ’ { c, [ -£ — 2- — 1 . 11 ] + c D j 


Cx,1 "«7”"' L ' “ c - (i i-i - *i> 


+ , U C tan “o ^ - >i-l) ; -1/2 

U c - (X,.! - *,) 


(9a) 


and 


U 2 . 

C y,i ■ 7T < C L + C D C 
U C 


j c tana o * - Vi ) 

U C - <*1-1 - *i> 


]} 


+ r u c tan °Q + - ?1-1> -1/2 


(9b) 


Use of Equation (5) into Equations (9a) and (9b) yields the desired 
expressions for the force coefficients as [9] 


c x,i - TT? < [ u c ta % + »i - >1-1» 2 + - <*1-1 - V* 1/2 


{C L [U c tana 0 + (yj - y^j)] + C Q [U c - (x^j - x.,)]} 


(10a) 


and 


C y,i * TV {CU c ta % + <>1 - Vl^ + < U c * <*1-1 - *1>] 2 > 


{C L CU c - { *i-l - V ] + C D [U c tana o + ( *i ' Vl )]} * 


2 , 1/2 


(10b) 
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The excitation forces may now be written as 


F x,i ■ \ <r u c ta "“o + ( *i - Vi^ 2 + [u c - (i i-i - J i )32 > 


{C L [U c tana 0 + (y. - y i _ 1 ) 3 + C p [U c - (x..^ - x^]} 


2,1/2 


(Ha) 


and 


F y,l 


■ \ pA {[U c tana 0 + (y. - y 1 _ 1 )3 2 + [U c - (x.^ - x.)] 2 } 


2,1/2 


(lib) 


(C L C u c - (x i _ 1 - x..)] + C Q [U c tana 0 + (y i - y^)]} . 

It may be noted that in addition to fluctuations in the excitation force due 
to convolute motion, there is also a fluctuation exressed quantitatively due 
to vortex shedding that was suggested in a qualitative way previously by 
others. As the vortex moves across the region of highest pressure on the 
convolute tip, the pressure is reduced. After the vortex has been shed the 
high pressure region is reestablished, thus producing an "on-off-on" type of 
fluctuation in the force which is already varying with convolute motion. 

Equations (11a) and (lib) are to be incorporated as excitation forces in 
the vibration model that represents the dynamic force balance of the bellows. 

The Lumped Parameter Vibration Model 

The development of an expression for the excitation force acting on an 
individual convolute suggests the need for a mechanical model which allows 
both longitudinal and transverse motion of the convolute. In order to develop 
such a model the bellows is divided into discrete mass elements such as the 
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one shown in Figure 9. By considering longitudinal and transverse motion of 
this mass element and the elements adjacent to it, a two dimensional vibration 
model is developed. This is illustrated in Figure 10, where Fq jX and Fg^y are 
damping forces which represent such effects as sliding between plys of the 
bellows and internal frictional damping. Viscous damping is not shown in the 
illustration, since the damping effect produced by the fluid surrounding the 
convolute is already included in the expression for the excitation. k x and ky 
are material spring constants, m m the material mass of one element and m f the 
fluid mass of one element. These components of the lumped model are further 
defined elsewhere in this report. 

From the vibration model illustrated in Figure 10 a system of equations 
describing the motion of the convolutions must be developed. This system must 
include all of the lumped parameters and effects described above. The system, 
when solved, must also provide displacement as a function of time for each of 
the convolute tips-. 

The development of the equations governing the system behavior results in 
two times the number of convolutions coupled, nonlinear, ordinary, 
differential equations. These equations are derived by performing a force 
balance on the discrete mass element shown in Figure 11(a). The force balance 
gives 


F = F . + F n sgn(x-) 
x x,i D,x 3 v v 


+ k x (x Hl ' 2x i + x f-l> ' <ra m + m F )x 1 


which upon rearranging becomes 


(m m + m F> x 1 - V S 9"<V + k x (2x i * x itl - x 1-l> * F x,i 


( 12 ) 
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Figure 9. Convolute Representing One Discrete Mass. 
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Figure 10. Vibration Model. 
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Figure 11a. Force Balance Along x-direction. 



Figure lib. Force Balance Along y-direction. 
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The force balance in the transverse direction also may be performed on the 
element shown in Figure 11(b). This gives 

EF y = F y,1 + F D,y + * 2 »i + ^1-1> ’ K + "f)»1 

Once again, rearrangement gives 

K + "f^l - Vy^'V + k y (2y 1 - *1+1 - *1-1 > = F y,i • (13) 

Upon substitution of the previously determined expressions for F x j and Fy j 
the system governing equations become- 

<%, + ra f> x 1 - Vx^V + 2 k x x 1 - k x x f+l - k x x 1-l 
• \ pA{[U c tann 0 + (J-, - y f _j)] 2 + [U c - i,)] 2 } 1/2 (14) 

{C L [U c tana 0 + (y. - y^)] + C D [U C - (x.^ - x.)]} 

and 

<m n, + m f>*i ' F D,y s9n(i ’l ) + Vi * k y*1+l ' Vi-1 

* j pA {[U c tano 0 + (y i - + [U c - (x^j - x,)] 2 } 1 '' 2 (15) 

{C L CU c - (x i ;l - X.)] + C D [U c tana 0 + (y i - y i-:l )]} , 

where each value of i designates a specific convolute. This produces the 
previously mentioned system of coupled, nonlinear, ordinary, differential 
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equations. 

The solution of this system is sought under the assumption that the 
bellows is rigidly attached at each end and that the entire system is 
initially quiescent. 

Components of the Lumped Parameter Model 

The Fluid Added Mass 

As a convolute is displaced it also displaces a certain amount of the 

fluid flowing through the bellows. This implies that the mass of this fluid 

must be included in the governing equations as an inertial effect; hence, the 
•• •• 

mf x.. and mf y.. terms. Each convolution has an instantaneous configuration 
ranging between the two extremes shown in Figure 12(a). This suggests that 
the fluid added mass varies with instantaneous longitudinal position much more 
than transverse position. The expression for the fluid added mass will 
therefore be a function of the longitudinal position of adjacent convolutes. 
With the geometry defined in Figure 12(b) an expression may be developed for 
the fluid mass. For example, the instantaneous cavity width and length of the 
straight section are given by 


« = « 0 + (X, - x 1 _ 1 ) 

(16) 

* ' *0 *7 (x f-l - x i> • 

(17) 


With this information the areas Aj, A£, and A 3 may be determined and swept 
around the bellows to yield the volume of the fluid added mass. Multiplying 
this volume by the fluid density, p , yields the expression for the fluid 
added mass, m^, as 
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m f • ’“V S < 8 o + (x i - x i-l» 2 + < s o + (x i - X i- 1>> [4 0 + \ 5 0 


+ t + ij— (x.. - x..^)] + ( —yr )C y («o + ( x -j ~ x -j_^)) + ^3^} • (18) 

Here, D m is the mean bellows diameter. Since this expression is multiplied by 
•• •• 

x.. an< i y i in the governing equations it introduces additional nonlinear 
effects into the system of governing equations. 

The Material Mass 

Each discrete mass in the vibration model is composed of the material 
mass of one convolution as well as the fluid added mass. Having determined an 
expression for the fluid mass added to one convolute in motion, an evaluation 
of the material mass which composes one convolute may also be obtained. 

A cross-section of one convolute is illustrated in Figure 13. The area 
of this cross-section is given by 

AREA = [(2ir + 4)a + 2h]t (19) 

The volume of material may then be determined by sweeping this area around the 
bellows mean diameter to yield 

VOLUME = TrD m [( 2tt - 4) a + 2h]t . (20) 

The material mass of one convolute, rr^ , is determined by multiplying this 
volume by the material density, p m . This gives 

m m = *D m P m [(2* - 4)a + 2h]t . (21) 
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Figure 13. Convolute Geometry 


to Calculate Material 
Mass. 



Figure 14a. Configuration to Calculate 
Spring Constant, k x - 



Figure 15a. Geometry to Calculate 
Convolute Tip Strain 
Along x-direction. 


Figure 14b. Configuration to 
Calculate Spring 
Constant, k . 





Figure 15b. Geometry to Calculate 
Convolute Tip Strain 
Along y-direction. 



The Longitudinal Spring Constant, k x 

The longitudinal spring force is the force tending to restore the 
convolute to its equilibrium position when it has been deformed by a 
longitudinal displacement. This force is a result of the elasticity of the 
bellows material . In the vibration model the spring force is represented as a 
longitudinal spring constant, k x , multiplied by the displacement of the 
convolute relative to adjacent convolutes. 

An expression for the spring constant can be determined by using the 
strain energy method along with the theorem of Castigliano. From the geometry 
illustrated in Figure 14(a) and with a force applied in the longitudinal 
direction the strain energy of this configuration is given by 


U = 


F x 2 « 


3ir 

T 


2 ) 


3 2 

a .ir a. it a , t , a £ 

wq t wq + tmj + mq wq 


, it a tt a 

+ t wq + ?wq 


+ ( 


4+ 2) 




(2 + tt) 


2 2 
a z A tt at 

wq + ? wq } 


( 22 ) 


Here, U is the strain energy for the configuration illustrated in Figure 
14(a), E is the bellows material modulus of elasticity, and G is the shear 
modulus. Ai, 1]^, are the cross-sectional area and moment of inertia for the 
section based on the bellows inside diameter. A 2 , I 2 and A 3 , I 3 are based on 
the bellows mean diameter and outside diameter, respectively. 

The theorem of Castigliano states that the derivative of the strain 
energy with respect to an applied force is the displacement of the point of 
application of the force in the direction of the force. Therefore, the 
longitudinal displacement of the convolute tip illustrated in Figure 14(a) may 
be determined as 
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* _ 9U _ or r/ 3 tT o\ 9 t TT a j_ TT d a 

6 x " " ZF x U “ Z) WT l + 7 Wj + 7 BG7~ + 


2 ... . 3 

^ a «, ^ tt a j.ira . f 3* 0 x a 

+ 7FTJ + 7 WJ^ + 7 + ' ~7 + Z ' 8E 1 3 


( 23 ) 


2 2 

, in , » a £ . n a£ , 

+ (2 + TT) -gtf- + ) • 


Since the longitudinal spring constant is the force required to produce a 
given displacement, the spring constant may be determined as 


k = 


l r , 3ir 0N a A ir a j_ ir a j_ £ 

7 {( T ■ Z) 7 7 7^ 8GA^ 


a 2 £ 


tt a tt a 
FTTJ T 7 WK^ + 7 mr 2 


+ Tfyr + -tt h f » - + tt mnr 1 + ( — ^ ^ 2 ) 


8EI. 


(24) 


2 2 

_i_ / a . \ a & . tt ait i *1 

+ ( 2 + tt) + > 


The Transverse Spring Constant, ky 

The transverse spring constant may be determined in a manner anlagous to 
the method used for determining the longitudinal spring constant. The strain 
energy for the configuration illustrated in Figure 14(b) is found to be 



3 2 

ira , tt a , ir a J _a£ J _£ 

7 7 WK^ 7 3g7[ wq + WK^ 


+ ( 




a + Z — i 

8EA 3 4 8GA 3 i 


(25) 


By the theorem of Catigliano the transverse displacement is given as 
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( 26 ) 


“3 O 

- _ 9U _oc r it a ir a ir a + a£ 

y " 9F y y 1 4 8EI X 4 8EA 1 4 8GAj 8EI 2 


+ ^ X f 9” /)> 3 + * 3 j. " 3 n 

+ + 1 "T 4) ■5TTJ + 7 + 7 } • 


The transverse spring constant may now be determined as 


+ £ 

+ WK 

III. 


F 3 ? 

y . 1 f it a x i a it a ^a^t 

? { ? WL[ + + 8EIJ 


9tt . x a . it a . tt a , -1 


J. I ^ a\ ° _L ° J. a , 

” 1 t 4) wrj + •? -goj + 7 wq } 


SOLUTIONS IN THE TIME AND FREQUENCY DOMAINS 


(27) 


Solution of the Governing Equations In the Time Domain 

The bellows vibrational response to flow excitation is obtained by 
solving the system of governing equations for the displacements. Since an 
analytical solution to this system is not known a numerical solution to be 
executed by computer must be chosen. In problems where the range of 
integration is considerable, the numerical method employed must be stable 
and/or relatively stable to obtain an accurate solution. Also, in a vibration 
problem where the motion is to be studied over a substantial number of 
oscillations, a method having a small per-step truncation error should be used 
in order to minimize the cumulative error. 

Upon consideration of these requirements, a "predictor - modifier - 
corrector" method called Hamming's Method was selected as the solution 
technique [10]. Hamming's Method is both stable and relatively stable and 
results in small per-step truncation errors. 

The iteration algorithm for a differential equation of the form 
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y' = f(y,t) (28a) 

may be stated as the predicted value of y at t + 1: 

P( W =y t-3 + ^ [2y 't - y 't-l + 2y 't-2 ] (28b) 

The modified value 

M ( y t+1^ = p ( y t+l^ " S “ c ^ y t^ ^ 28c ^ 

The corrected value 

C <W ’ I {9y t - y t-2 + 3h ["'( y ' t+ l ) + 2y 't - (28d) 

The final value 

F ^ y t+1^ = C ^ y t+1^ + T2T ^ y t+l^ ' C ( y t+1^* ^ 28e ^ 


Here, h is the time step used in the iteration procedure. 

As can be seen from these equations Hamming's Method also has the 
advantage of being relatively efficient in terms of computing time, since 
there is only one function evaluation per time step. A disadvantage, however, 
is that Hamming's Method is not sel f- starting. It requires solution values at 
the first three time steps in order to begin the iteration process. 

To provide the necessary solution for the first three time steps, a 
fourth-order Runge-Kutta procedure is employed. While being very inefficient 
in terms of computing time, the Runge-Kutta procedure does have the advantage 
of being self-starting. 
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By starting the iteration with an initially quiescent system, the Runge- 
Kutta procedure determines the longitudinal and transverse displacements of 
each convolute tip for the first three time steps. This information is then 
passed on to Hamming's Method which iterates to 4096 time steps, each step 
being of 0.01 millisecond duration. The displacement information for the last 
2048 time steps is stored as a "data sample" for frequency and fatigue 
calculations. By choosing the last 2048 time steps the "data sample" is 
considered far enough from time zero where steady-state vibration is assumed. 

With the 0.02048 second sample of steady-state vibration displacement 
information, a frequency response analysis and strain analysis may be 
performed. Plots of displacement versus time from this solution are shown in 
Appendix 1. 

Frequency Response 

Examination of the displacement versus time plots reveals that both 
longitudinal and transverse vibrating of an individual convolute are made up 
of several frequency components. Each of these frequency components 

contribute to the damage done to the convolute during a given time period. In 
order to estimate the extent of this damage on each convolute the dominating 
frequencies of vibration in both the longitudinal and transverse directions 
must be determined. 

With the time-dependent displacement obtained from the solution of the 
governing equations, a Fourier Analysis may be performed. This analysis 
yields the frequency response of each convolute in each coordinate 

direction. The Fourier Analysis is performed by means of the Fast Fourier 
Trans from algorithm which is developed here. 

By representing the convolute response as a complex Fourier series one 
obtains 
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x(t) = 


oo 


(29a) 



(29b) 

(30a) 



(30b) 


to n is a discrete angular frequency and t is time. Recall from the solution of 
the governing equations that the sampling period, T, and the number of data 
samples, N, are respectively. 


T = 0.02048 second 


N = 2048 samples 


Therefore the time step size is 


Ola) 



The fundamental frequency, oj q , is given by 



(31b) 
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(31c) 



and the maximum frequency which may abe determined, u> 


_ 2ir _ 2ir _ M 
“max " At ' TTTnT " N “o * 


(33d) 


Define a discrete set of N frequencies at which the convolute response 
will be determined. 


a) = n«o = — j-; n = 0, 1, 2, ...» N 


Since x(t) and y(t) vanish for time less than zero and since they have not 
been determined for time greater than T, the transform may be approximated as 


X(o>) = f x(t)e'^ a>t: dt 

J o 

Y(u>) = J y(t)e' 1u)t dt 


(33a) 


(33b) 


If these integrals are approximated by summations, the expressions become 


II- JL 

E «< v 


(34a) 


») -£y(V 


(34b) 


Upon evluating these transform expressions at discrete frequencies, io n , one 


obtains 


X(»„> 


l E x(t m exp[ ' i T 21 ] 

m=0 


( 35a) 
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and 


N-l 

V <“r>> ’ i * 

m=0 

Let = exp[-i ] to get 

N-l 

and m=0 

N-l 

V K> -l E ''W™ • 

m=0 

Each power of increases the polar angle by 2ir/N . After N 
pattern repeats so it is not necessary to compute all powers, 
define such that 

U k = (W N ) k ; 0 < k < N 
and 

U k = U mod(k,N) * k > N » 

where M0D(k,N) is the modulus function. The combination of 
previous transform expression yields the equations for the 
response at discrete frequencies as 

N-l 

X (“n) = IT U M0D(mn,N) 

and m=0 

N-l 

Y ^“n^ = 7T / a U M0D(mn,N) * 

m=0 


(35b) 


(36a) 

(36b) 

powers the 
Therefore, 


(37) 

with the 
transformed 

(38a) 

(38b) 
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Note that since (W n ) m ^“ n ^ is the complex conjugate of (Wn|) n,n , it is sufficient 
to evaluate X(<D n ) and Y(a> n ) for n > 0 . 

The results of the Fast Fourier Transform are the frequency response 

plots shown in Appendix 1. With this information the dominant frequencies of 
both longitudinal and transverse vibration of each convolute may be 
determined. These dominant frequencies will later be used to assess the 
cumulative damage effects on each convolute. 

IV. CYCLIC STRAIN, DAMAGE AND FATIGUE LIFE 
Estimation of Cyclic Strains 

The numerical model to calculate fatigue life of bellows is based on 

obtaining strain amplitude for the cyclic deformation of each convolute. The 
strain amplitude has to be evaluated at the convolute tip since this is the 
location at which the bellows is known to fail. A maximum and minimum value 
for the tip strain resulting from both longitudinal and transverse vibration 
needs to be determined. Since the strain induced at a point by different 
mechanisms is additive (same fibres damaged), the maximum tip strain resulting 
from longitudinal vibration and the maximum tip strain resulting from 

transverse vibration may be added to get a total maximum convolute tip 
strain. Likewise, the corresponding minimum strains may be added to get a 
total minimum convolute tip strain. The strain amplitude may then be 

determined by the equation: 


Ae 


e max “ e mi n 
2 


where e is the tip strain and Ae the strain amplitude. 


( 39 ) 
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In order to determine an expression for the convolute tip strains 
resulting from both longitudinal and transverse displacements the strain 
energy method along with the theorem of Castigliano is used. For longitudinal 
displacement of the convolute tip it is assumed that the strain energy 
resulting from the displacement is concentrated in the section of the 
convolute illustrated in Figure 15(a). The strain energy for this 
configuration is given by 


■ <ViV 2 « 


3tt 

“5 


1 ) 




"5 T^T 




1 ta 3 , 1 ta , 
7TI^ + 71^ } ’ 


(40) 


where, as before, E is the bellows material modulus of elasticity and G is the 
shear modulus. Aj and Ij are the cross-sectional area and the moment of 
inertia based on the bellows inside diameter. A 2 , 1 2 and A 3 , I 3 are based on 
the bellows mean diameter and outside diameter, respectively. 

By applying the theorem of Castigliano the longitudinal displacement of 
the tip, x-j , may be determined as 


x i 


au 


3ir 




tt a , 1 ta , lia , 

8 GA X 2 EI 2 2 GA 2 1 * 


Since a . = Ee . , the displacement may also be expressed 

X j i X j 1 


(41) 


3 7T 




( 42 ) 
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1 ta 3 . 1 ta 3 1 ta . 

7^ + 7HJ + 7^“ } 


and the convolute tip strain may be written as 


x,i 


x. ~ ,3 , 

1 r /3ir . x a , it a , ir a 

WK[ { [ S' TT^ + S1^E + SW[ 


, 1 ta 3 . 1 ta 3 . 1. ta 
+ ?^ + ?TI^ 7 J GA^ 


(43) 


-1 


Since the transverse displacements are smaller than the longitudinal 
displacements the strain energy resulting from transverse displacements is 
considered to be concentrated in a smaller section of the convolute. Figure 
15(b) is an illustration of this section. Once again the strain energy for 
this configuration is given by 


* <ViV 


TT 




T3I 


TSSJ' 


(44) 


and application of Castigliano' s theorem gives 


y i = 


20 


7KT = 2(a y»i A i ) 


y,i l 


* a 3 a a 

* { + + wq 


(45) 


Upon substituting a . = e .E and solving the resulting expression for e ^ , 
the equation for the convolute tip strain resulting from transverse vibration 
is obtained as 


_ y i 4 , a 3 a a -1 
e y,i ~ { TT[ + TK[ + 


(46) 
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As described earlier, the maximum values of e Y * and e„ _• which occur during 

the sampling period are added to get a total maximum tip strain, e. _ . The 

i ,max 

minimum values of e . and e . are also added to yield a total minimum tip 

X , 1 jj ‘ 

strain, e- ■ . The convolute tip strain amplitude, Ae^ , may then be 

I jHI III I 

determined by 


Ae i = 7 


( C • 

v i ,max 


e i ,min^ 


(47) 


Estimate of Cycles to Failure 

A modern approach to characterize the fatigue behavior of materials is to 
focus on the cyclic strain life [12], The total strain, e , is considered as 
having an elastic e Q , and a plastic, Ep , components. Expressed as strain 
amplitudes, this implies 


As 

T 




(48) 


The elastic strain-life can be expressed as 

A s - 0 ^ i u 

“F " *r " r a f { 2 N f } * (49) 

i 

where and , respectively, represent the true stress amplitude and the 
fatigue strength coefficient, E is the modulus of elasticity, b is the fatigue 
strength exponent (typically, b varies between -0.05 and -0.12) and is the 
number of cycles to failure. The plastic strain-life is related by a power 
function as 



( 50 ) 
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I 

where is the fatigue ductility coefficient and d is the fatigue ductility 
exponent (typically, d varies between -0.5 and -0.7). The total strain 

amplitude may thus be expressed as 

o 1 

(2N t ) b + e f ' (2N f ) d . (51) 

This equation is called the strain-life relation and forms the basis for the 
strain-life approach to predicting fatigue behavior of such material as 
wrought metals. For thin-walled shells, such as the bellows, Nf is typically 
taken to be the number of cycles to initiate a crack. 

With the strain amplitude for each convolute in hand it is possible to 
compute the number of cycles to failure for each convolute. For the types of 
steels from which bellows are commonly made, the equation constants are given 
by 


b - -0.1 


0 ^' « 150 kpsi 

0.5 < e f ' < 1 (52) 

and d « -0.7 

The strain-life relation obviously requires an iterative solution. To 
accomplish this, simple Newton-Raphson Method has been employed. Since the 
slope of the strain-life relation is always negative, there should be no 
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problems with convergence so long as N f is initially relatively small in the 
iteration process. 

Cumulative Damage and Fatigue Life of Bellows 

By estimating the cummulative damage done to each convolute during the 
sample period an estimation can be made as to which convolute will fail 
first. The assumption made here is that the convolute undergoes cyclic strain 
of amplitude Ae.. for the time of the sampling period at each of the dominant 
frequencies determined by the Fourier Analysis. For example: 


o convolute #3 vibrates with strain amplitude Ae^ at a frequency fj 
for a time equal to the sampling period, 
o convolute #3 then vibrates with strain amplitude Aeg at frequency f 2 
for a time equal to the sampling period 
o etc. 

From this approach the cumulative damage done to each convolute during 
one sampling period can be determined. This is done by use of the Palmgren- 
Miner linear-cumulative-damage rule [12] which may be stated as 


2Nj (Reversals applied at a al ) 
d " ZN f ^ " (Reversals to failure at a al ) * 

where d = damage. Failure will occur when 


(53) 



(54) 


In the present model the damage for each convolute during the sampling 
period may be determined by 
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+ H 

Here T is the sampling period, Nfj is the cycles to failure for the i 
convolute, f ^ . f ^ n are the dominant frequencies of vibration, and d^ is the 
damage for the i^ convolute. Since steady state vibration is assumed to 
exist, the largest value of d^ corresponds to the convolute which is likely to 
fail first. 

The number of cycles to failure for this convolute is also the cycles to 
failure for the entire bellows. A failure time window may be determined by 
dividing the bellows cycles to failure by the maximum and minimum dominant 
frequencies of the convolute which is first to fail. This defines a minimum 
time to failure as 


Snin “ f. mav 
i,max 


( 56 ) 


and a maximum time to failure as 


t = Nf 
max Ti~7n * 


(57) 


This t min and t max then establish a failure time window. 



V. CLOSING REMARKS 


In testing the bellows computer program it became apparent that the best 
results are obtained for bellows with up to five convolutions. For bellows 
with greater than five convolutions the cycles to failure begins to decrease 
considerably. Therefore, it is recommended that program BELLOWS be used for 
bellows with up to five convolutions. 

The limitation of the program to give resonable vibration amplitudes for 
bellows with less than five convolutions is attributed to cumulative 
displacement effects in the springs connecting the discrete masses of the 
mechanical model. For more than five convolutions the displacements in each 
consecutive spring adds to eventually produce unreasonable convolute tip 
displacement. Rather than introducing artificial damping into the governing 
equations to control these amplitudes, it was decided to limit the model to a 
five convolution model. Additional effects due to stress wave propagations in 
bellows with large number of convolutions need to be examined further in 
future research. 

For bellows with differing number of convolutions the interaction between 
the convolutes results in the strain energy being concentrated in a smaller or 
larger volume of the convolute. Fewer convolutions result in the strain 

energy being concentrated in a smaller volume near the convolute tip. To 
incorporate this effect into the model a factor multiplying the equations for 
strain calculation has been introduced. The strain multiplying factor varies 
linearly with the number of convolutions and its determination is internal to 
the computer program. 

A thorough literature search into the effects of sliding damping between 
plys and internal frictional damping did not yield sufficient information to 
explicitly incorporate these effects into the model. Just the same such 
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effects do exist and were included in the model via damping forces. A 
comparison of calculated fatigue life cycles with measured values [2] indicate 
very encouraging agreements. The computer calculations were made for 
equil valent test conditions and bellows specifications to those used in the 
experiments, except that the number of convolutions used in the model was less 
than or equal to five. The agreement between the two indicated in Tables 1 
and 2 is within an order of magnitude. It is concluded therefore that within 
the scope of its application, the model presented in this study is quite 
sati sfactory. 


Table 1 Computer Results 


Bellows 

Number 

# Convolutes 
in Model 

Flow 

Velocity (ft/s) 

Cycles 

to Fail (10 6 ) 

Min. Time 
to Fail 
( sec) 

Max. Time 
to Fail 
( sec) 





5028 

5 

35 

2.9 

979. 

11,747. 

5028 

5 

50 

210. 

77,851. 

856,362. 

5028 

5 

65 

6.7 

2482. 

27,301. 

5013 

3 

35 

0.67 

651. 

4556. 

5013 

3 

50 

0.056 

63.8 

383. 

5005 

3 

35 

0.33 

285. 

1708. 

5005 

3 

50 

0.025 

24.5 

171.2 


Table 2. MSFC Tests [2] 


Bellows 

Flow 

Cycles to 

Time to 

Number 

Velocity (ft/s) 

Fail (10 6 ) 

Fail (sec) 





5028 

65.8 


787. 

5028 

66.6 

.16 

150. 

5013 

50 

1.3 

2007 - 2027 

5013 

45 

.47 - .48 

749 - 760 

5005 

33.2 

.78 - .79 

1841 - 1863 

5005 

65.2 


75 - 100 
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APPENDIX I 


SAMPLE PLOTS FOR TIME AND FREQUENCY DOMAIN RESPONSES 
In the time response plots, the core velocity, U c , is taken as 45 
ft/ sec. The horizontal axis indicates time in seconds and the vertical axis 
indicates displacement in inch, x and y, respectively, indicate longitudinal 
and transverse vibrations displacements. Subscripts on x and y indicate the 
data for the particular convolution number. 

In the frequency response plots, the vertical axis indicates the Fourier- 
transformed displacements and the horizontal axis incidates frequency in 
radians/ second. The core velocity, U c , for all plots is 45 ft/ sec. 
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APPENDIX II 


LISTING OF PROGRAM BELLOWS 
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o o n o non 


PROGRAM BELLOWS (INPUT , OUTPUT , TAPE5=INPUT , TAPE6-0UTPUT) 
C 

C ***************************************************** 


C * * 

C * BELLOWS FATIGUE PROGRAM * 

C * 1/28/86 * 

C * * 

C * * 

C * * 

C * INVESTIGATOR: DR. P. V. DESAI * 

C * ASSISTANT: L. D. THORNHILL * 

C * * 


Q ***************************************************** 

C 

C THIS PROGRAM COMPUTES AN ESTIMATE OF THE FATIGUE 
C LIFE OF BELLOWS UNDER A VARIETY OF OPERATIONAL 
C CONDITIONS. 

C 

CALL ENTER 
CALL RUNGE 
CALL HAMM 
CALL STRAIN 
CALL FFT 
CALL CYCLE 
CALL DAMAGE 
CALL RESULT 
C 

STOP 

END 


SUBROUTINE ENTER 


ENTER GEOMETRIC, MATERIAL AND FLOW PARAMETERS AND 
CALCULATE EQUATION CONSTANTS. 


C 


COMMON/BLOK1 /PI, DELO , P IDEN , ALO , TE , EM , UC , UCT , N 
COMMON / B LOK2 / XK , YK , FDX , FDY , DENSA 
COMMON/BLOK5/ Al , E , BOTX2 , BOTY2 
DIMENSION FDX (9) , FDY (9) 


PRINT*, 'ENTER 
READ*, PITCH 
PRINT*, 'ENTER 
READ*, HI ' 
PRINT*, 'ENTER 
READ*, TPLY 
PRINT*, 'ENTER 
READ*, NPLY 
PRINT*, 'ENTER 
READ*, DIN 
PRINT*, 'ENTER 
READ*, N 
PRINT*, 'ENTER 
READ*, DENF 
PRINT*, 'ENTER 
READ* , UFEET 
PRINT*, 'ENTER 
READ*, El 


CONVOLUTE PITCH IN INCHES.' 
CONVOLUTE HEIGHT IN INCHES.' 

PLY THICKNESS IN INCHES. ' 

NUMBER OF PLYS. ' 

INSIDE BELLOWS DIA. IN INCHES.' 
NUMBER OF CONVOLUTES.' 

FLUID DENSITY IN LBM/ (IN**3) . ' 
FLOW VELOCITY IN FT/S.' 

THE ELASTIC MODULUS IN PSI. ' 
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</>{/></> 


PRINT*, 'ENTER THE MATERIAL DENSITY IN LBM/ (lN**3) . ' 
READ*, DENM 
C 

TE=TPLY*FLOAT (NPLY) 

UOUFEET*12. 

E=E1*12. 

G-E/ (2.*(l.+l./3.)) 

C 

PRINT*, 'CONVOLUTE PITCH IS '.PITCH,' IN.' 

PRINT*, 'CONVOLUTE HEIGHT IS ',HI,' IN.' 

PRINT*, 'PLY THICKNESS IS ' , TPLY, ' IN.' 

PRINT*, 'NUMBER OF PLYS IS ' ,NPLY 

PRINT*, 'INSIDE BELLOWS DI A. IS ' ,DIN, ' IN.' 

PRINT*, 'NUMBER OF CONVOLUTES IS ' ,N 

PRINT*, 'FLUID DENSITY IS ' ,DENF, ' LBM/(IN**3)' 

PRINT*, 'FLOW VELOCITY IS ' ,UFEET, ' FT/S' 

PRINT*, 'ELASTIC MODULUS IS ',E1,' PSI' 

PRINT*, 'MATERIAL DENSITY IS ' ,DENM, ' LBM/(IN**3)' 

C 

DO 3 J-1,N 
FDX(J)-8. 

FDY (J)=2. 

3 CONTINUE 
C 

C CALCULATE EQUATION CONSTANTS 
C 

IF <UC .GT. 0. .AND. UC .LE. 180.) ALPHA- 1.3 
IF (UC .GT. 180. .AND. UC .LE. 300.) ALPHA=1.2 

IF (UC .GT. 300. .AND. UC .LE. 420.) ALPHA-1. 0 

IF (UC .GT. 420. .AND. UC .LE. 540.) ALPHA-0.6 

IF (UC .GT. 540. .AND. UC .LE. 660.) ALPHA-0.4 

IF (UC .GT. 660. .AND. UC .LE. 780.) ALPHA-0.1 

IF (UC .GT. 780.) ALPHA-0.05 
R- PITCH/4. 

PI- 4.*ATAN(1.) 

DOUT-DIN+2 . *HI+2 . *TE 
DM- (DIN+DOUT) /2. 

SIGMA- (2 . *R) +TE 
DELO- (2 . *R) -TE 
ALO- HI- (2. *R) 

PIDEN- PI*DM*DENF 

EM- PI*DM*DENM* (R* ( (2 . *PI) -4 . ) +2 . *HI) *TE 
C 

AMOM1- (PI*DIN*(TE**3.))/12. 

AMOM2- (PI*DM*(TE**3.))/12. 

AMOM3- (P I *DOUT* (TE**3 . ) ) / 1 2 . 

C 

Al=PI*DIN*TE 

A2=PI*DM*TE 

A3-PI*DOUT*TE 

B=HI-2.*R 

C 

BOTTOM- (PI/32 . ) * ( (R**3 . ) / (E*AM0M1) + R/(E*Al) 

+ R/(G*A1) + R/(E*A3) + R/(G*A3)) + 0. 125* (( (R**2.) / 
(E*AMOM2))*B + B/(E*A2)) + ((9.*PI-16.)/32.)*((R**3.)/ 
(E*AMOM3) ) 

B0TY2- (PI/32 . ) * ( (R**3 . ) / (E*AM0M1) +R/ (E*Al) +R/ (G*A1) ) 

C 

YK- 1. /BOTTOM 
C 
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Cl=(3.*PI)/8.-l. 

C2=PI/8. 

C3=1.+(3.*PI) /8. 

C4=l .+PI/2. 

C5=Pl/4. 

C 

B0TX“C1* (R**3 . ) / (E*AMOMl) + (C2*R* (A1+A3) * (E+G) ) / (Al*A3*E*G) 
$ + ( (R**2 . ) *B) / (2 . *E*AM0M2) +B/ (2 . *G*A2) +C3* (R**3 .) / (E 

$ *AM0M3) + (C4* (R**2 . ) *B) / (E*AM0M3) + (C5*R* (B**2 . ) ) / (E 

$ *AM0M3) 

BOTX2=BOTX-C3* (R**3 . ) / (E*AM0M3) -C4* (R**2 . ) *B/ (E*AM0M3) 

$ -C5*R* (B**2 . ) / (E*AM0M3) -C2*R/ (E*A3) -C2*R/ (G*A3) 

C 

XK“1./B0TX 

A= (R+TE/2 . ) * (P 1**2 . ) *DIN 
DENSA= 0.5*DENF*A 
UCT= UC*TAN (ALPHA) 

C 

RETURN 

END 

C 

C 

C 

SUBROUTINE RUNGE 
C 

C CALCULATE STARTING VALUES FOR HAMMINGS METHOD. 

C 

COMMON/BLOK1/ PI,DELO,PIDEN,ALO,TE,EM,UC,UCT,N 
COMMON/BLOK3/ X,Y 

COMMON/BLOK35/ XP , XP 1 , XP2 , XP3 , XP4 , YP , YP1 , YP2 , YP3 , YP4 , 

$ XPP , XPP 1 , XPP2 , XPP3 , XPP4 , YPP , YPP 1 , YPP2 , YPP3 , YPP4 , 

$ CD, CL 

DIMENSION X(0:6, 0:2048) ,Y(0:6, 0:2048) ,XP(0:6) ,XP1(0:6) , 

$ XP2 (0:6) ,XP3(0:6) ,XP4(0:6) ,YP(0:6) ,YP1(0:6) , 

$ YP2 (0:6) ,YP3(0:6) ,YP4(0:6) , XPP (0:6) ,XPP1(0:6) , 

$ XPP2 (0:6) ,XPP3 (0:6) ,XPP4(0:6) , YPP (0:6) ,YPP1 (0:6) , 

$ YPP 2 (0:6), YPP3 (0 : 6) , YPP4 (0 : 6) , CD (5) , CL (5) 

COMMON / B LOK4/ TX , TXP , TY , TYP , K, Q , H , L , J , I 
REAL K 

DIMENSION TX (3) , TXP (3) , TY (3) , TYP (3) , K (4) , Q (4) 

C 

C TIME STEP, H, AND CONSTANTS. 

C 

H= 0.00001 
Al- H/6. 

A2= 1./6. 

A3= H/2. 

A4= H/4. 

C 

C INITIAL CONDITIONS 
C 

DO 10, J“0,N+1 
X(J,0)=0.0 
Y(J,0)=0.0 
XP(J)=0.0 
XP1(J)=0.0 
XP2(J)=0.0 
XP3(J)=0.0 
XP4(J)=0.0 
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YP(J)=0.0 

YPl (J) "0.0 

YP2(J)=0.0 

YP3(J)=0.0 

YP4(J)=0.0 

XPP1(J)=0.0 

XPP2(J)=0.0 

XPP3(J)=0.0 

XPP4(J)=0.0 

YPP1 (J) “0.0 

YPP2(J)=0.0 

YPP3(J)=0.0 

YPP4(J)=0.0 

10 CONTINUE 
C 

C BOUNDARY CONDITIONS 
C 

DO 30, J=0,N+1,N+1 

DO 20, 1-0,2048 
X(J,l)=0.0 
Y (J, I) =0.0 

20 CONTINUE 

30 CONTINUE 
C 

C RUNGE KUTTA STARTER. 

C 

DO 50, 1=0,3 
C 

DO 31, J=1,N 
XP4 ( J) =XP3 (J) 
XP3(J) =XP2(J) 

XP2 (J) =XP1(J) 
XP1(J)=XP(J) 
YP4(J)=YP3 (J) 
YP3(J)=YP2 (J) 

YP2 (J)=YPl (J) 

YPl (J)=YP (J) 

XPP4 ( J) =XPP3 ( J) 
XPP3 ( J) =XPP2 ( J) 
XPP2 ( J) =XPPl ( J) 
YPP4 ( J) =YPP3 ( J) 
YPP3(J)=YPP2(J) 
YPP2 (J) =YPPl (J) 

31 CONTINUE 
C 

DO 40, J=1,N 
C 

DO 32, M=l,3 
IND= J+(2-M) 

TX (M) =X (IND, I) 
TXP (M) =XP1 (IND) 
TY (M) =Y (IND, I) 
TYP (M) =YP1 (IND) 

32 CONTINUE 
C 

L=1 

CALL KQ 
C 

DO 34, M=l,3 
IND=J+ (2-M) 



TX (M) =X (IND , I) +A3*XP 1 (IND) 

TXP (M) =XP 1 (IND) + (K ( 1) /2 . ) 

TY (M) =Y (IND, I) +A3*YP1 (IND) 

TYP (M) =YPl (IND) + (Q (1) /2 . ) 

34 CONTINUE 
C 

L=2 

CALL KQ 
C 

DO 36, M-=l,3 
IND=J+(2-M) 

TX (M) =X (IND, I) +A3*XP1 (IND) +A4*K ( 1) 

TXP (M) =XP 1 (IND) + (K (2) /2 . ) 

TY (M) =Y (IND, I) +A3*YP1 (IND) +A4*Q (1) 

TYP (M) =YP1 (IND) + (Q (2) /2 . ) 

36 CONTINUE 
C 

L=3 

CALL KQ 
C 

DO 38, M=1 , 3 
IND=J+ (2-M) 

TX (M) =X (IND, I) +H*XP 1 (IND) +A3*K(2) 

TXP (M) -XP 1 (IND) +K (3) 

TY (M) -Y (IND, I) +H*YP 1 (IND) +A3*Q (2) 

TYP (M) =YP1 (IND) +Q (3) 

38 CONTINUE 
C 

L=4 

CALL KQ 
C 

X ( J , 1+1) =X ( J , I) +H*XP1 ( J) +A1* (K (1) +K (2) +K (3) ) 

XP ( J) -XP 1 ( J) +A2* (K ( 1) +2 . *K (2) +2 . *K (3) +K (4) ) 

Y ( J , 1+1) =Y ( J , I) +H*YP 1 ( J) +A1* (Q (1) +Q (2) +Q (3) ) 

YP ( J) =YP 1 ( J) +A2* (Q (1) +2 . *Q (2) +2 . *Q (3) +Q (4) ) 

C 

40 CONTINUE 
50 CONTINUE 
C 

RETURN 

END 

C 

C 

C 

SUBROUTINE KQ 
C 

COMMON/BLOK1 /PI, DELO , PIDEN , ALO , TE , EM, UC , UCT , N 
COMMON/BLOK2/ XK , YK , FDX , FDY , DENSA 
DIMENSION FDX (9) , FDY (9) 

COMMON/BLOK3/ X,Y 

COMMON/BLOK35/ XP , XP 1 , XP2 , XP3 , XP4 , YP , YP 1 , YP2 , YP3 , YP4 , 

$ XPP , XPP 1 , XPP2 , XPP3 , XPP4 , YPP , YPP 1 , YPP2 , YPP3, YPP4 , 

$ CD, CL 

DIMENSION X(0:6, 0:2048) ,Y(0:6, 0:2048) ,XP(0:6) ,XP1(0:6) , 
$ XP2 (0:6) ,XP3(0:6) ,XP4(0:6) ,YP(0:6) ,YP1 (0:6) , 

$ YP2 (0:6) , YP3(0:6) ,YP4(0:6) ,XPP(0:6) ,XPP1(0:6) , 

$ XPP2 (0:6) ,XPP3 (0:6) ,XPP4(0:6) ,YPP(0:6) ,YPP1(0:6) , 

$ YPP2 (0:6), YPP3 (0:6), YPP4 (0:6), CD (5) , CL (5) 

C0MM0N/BL0K4/ TX , TXP , TY , TYP , K , Q , H , L , J , I 
REAL K 
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DIMENSION TX(3) ,TXP(3) ,TY(3) ,TYP(3) ,K(4) ,Q(4) 

C 

DI FF-DELO+TX (2) -TX (3) 

XMAS- PIDEN* ( (PI/8 . ) * (DIFF**2 . ) +DIFF* (ALO+O . 5*DEL0+TE 
$ +0 . 5* ( 1 . -PI) * (TX (2) -TX (3) ) ) +0 . 5* (4 . -PI) * ( (0 . 5*DIFF 

$ +TE)**2.))+EM 

YMAS-XMAS 
C 

CD(1)— 0.5 
CL(l)-0.2 
DO 10 KOUNT-2 , N 
CD (KOUNT) -0.085 
CL (KOUNT) “0 . 08 
10 CONTINUE 
C 

TERMl-UCT+TYP (2) -TYP (3) 

TERM2-UC+TXP (2) -TXP (3) 

C 

XPP1 (J) =FCN1 (XMAS , TERM1 , TERM2 , TX (1) , TX (2) , TX (3) , TXP (1) , 
$ TXP (2) , TXP (3) , TYP (2) , TYP (3) , CD ( J) , CL ( J) , FDX ( J) , FDY ( J) ) 
K(L)=H*XPP1(J) 

C 

YPPl (J)“FCN2 (YMAS, TERMl, TERM2,TY(1) ,TY (2) ,TY(3) ,TYP(1) , 
$ TYP (2) , TYP (3) , TXP (2) , TXP (3) , CD ( J) , CL ( J) , FDX ( J) , FDY ( J) ) 
Q(L) “H*YPP1 (j) 

c 

RETURN 

END 


FUNCTION FCNl (XMAS , TERM1 , TERM2 , Xl , X2 ,X3 , VXl , VX2, VX3 , VY2 , 
$ VY3 , CD , CL , DX , DY) 

COMMON / B LOK2 / XK , YR , FDX , FDY , DENS A 
DIMENSION FDX (9) , FDY (9) 

FCN1“ (DENSA* (SQRT ( (TERMl**2 . ) + (TERM2**2 . ) ) ) * (CL*TERMl 
$ +CD*TERM2) +DX*SIGN (1 . , VX2) +XK* (Xl-2 . *X2+X3) ) /XMAS 

RETURN 

END 


FUNCTION FCN2 (YMAS , TERMl , TERM2 , Y1 , Y2 , Y3 , VYl , VY2 , VY3 , VX2 , 
$ VX3,CD,CL,DX,DY) 

COMMON / B LOK2 / XK , YR , FDX , FDY , DENSA 
DIMENSION FDX (9) , FDY (9) 

FCN2- (DENSA* (SQRT ( (TERMl **2 . ) + (TERM2**2 . ) ) ) * (CL*TERM2 
$ +CD*TERM1) +DY*SIGN (1 . , VY2) +YR* (Yl-2 . *Y2+Y3) ) /YMAS 

RETURN 

END 


SUBROUTINE HAMM 

IMPLEMENTING HAMMINGS METHOD FOR SOLVING GOVERNING EQUATIONS. 
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COMMON/B LOK1 / P I , DELO , P IDEN , ALO , TE , EM , UC , UCT , N 
COMMON/BLOK2/ XK , YK , FDX , FDY , DENSA 
DIMENSION FDX (9) , FDY (9) 

COMMON / BLOK3 / X,Y 

COMMON/BLOK35/ XP , XP 1 , XP2 , XP3 , XP4 , YP , YPl , YP2 , YP3 , YP4 , 

$ XPP , XPP 1 , XPP2 , XPP3 , XPP4 , YPP , YPP 1 , YPP2 , YPP3, YPP4 , 

$ CD CL 

DIMENSION X(0:6, 0:2048) ,Y(0:6, 0:2048) ,XP(0:6) ,XP1(0:6) , 
$ XP2 (0 : 6) , XP3 (0 : 6) , XP4 (0 : 6) , YP (0 : 6) , YP 1 (0 : 6) , 

$ YP2(0:6) ,YP3(0:6) ,YP4(0:6) ,XPP(0r6) ,XPP1(0:6) , 

$ XPP2 (0:6) ,XPP3(0:6) ,XPP4(0;6) , YPP (0:6) , YPP 1(0: 6) , 

$ YPP2 (0 : 6) , YPP3 (0 : 6) , YPP4 (0 : 6) , CD (5) , CL (5) 

DIMENSION PX (0:6) ,PX1(0:6) ,PXP(0:6) ,PXP 1(0:6) ,PY(0:6) , 

$ PY1 (0:6) ,PYP(0:6) ,PYP1(0:6) ,CX(0:6) , CXI (0:6) ,CXP(0:6) , 

$ CXP1 (0:6) ,CY (0:6) ,CY1 (0:6) ,CYP(0:6) ,CYP1 (0:6) ,MX(0:6) , 

$ MXP(0:6) , MY (0:6) ,MYP(0:6) 

REAL MX,MXP,MY,MYP 
C 

C TIME STEP, H, AND EQUATION CONSTANTS. 

C 

H=0. 00001 
A«(4.*H)/3. 

B=112./121. 

C=l./8. 

D=3.*H 

AB=9./121. 

C 

C INITIAL VALUE LOOP. 

C 

DO 3, J=0,N+1 
1=3 

PX(J) =X(J,I) 

CX(J)=X(J,I) 

PXP(J)=XP1(J) 

CXP(J)=XP1(J) 

PY (J)=Y (J, I) 

CY (J) =Y (J, I) 

PYP(J) =YP1 (J) 

CYP(J)=YPl (J) 

3 CONTINUE 

C 

C BOUNDARY VALUE LOOP. 

C 

DO 9, J=0,N+1,N+1 
MX (J) =0.0 
MXP (J)=0.0 
MY (J) =0.0 
MYP (J)=0.0 
9 CONTINUE 

C 

C *** HAMMINGS METHOD *** 

C 

DO 60, K=1 ,2 
C 

C REESTABLISH INITIAL CONDITIONS. 

C 

IF (K .EQ. 1) GO TO 14 
DO 12, 1=0,3 

IND=2048-3+I 7 , 



DO 10, J=0,N+1 
X ( J , I) =X ( J , IND) 

Y (J, I) =Y (J, IND) 

10 CONTINUE 
12 CONTINUE 
C 

14 DO 50, 1=3,2047 
DO 20, J=1,N 
C 

C PREDICTOR 
C 

PXP 1 ( J) =XP4 ( J) +A* (2 . *X?P 1 ( J> -XPP2 ( J) +2 . *XPP3 ( J) ) 

PX1 ( J) =X ( J , 1-3) +A* (2 . *XP1 ( J) -XP2 ( J) +2 . *XP3 ( J) ) 

PYP 1 ( J) =YP4 ( J) +A* (2 . *YPP1 ( J) -YPP2 ( J) +2 . *YPP3 ( J) ) 

PY1 ( J) =Y ( J , 1-3) +A* (2 . *YP 1 ( J) -YP2 ( J) +2 . *YP3 ( J) ) 

C 

C MODIFIER 
C 

MXP (J) =PXP1 ( J) -B* (PXP ( J) -CXP ( J) ) 

MX ( J) =PX1 ( J) -B* (PX ( J) -CX ( J) ) 

MYP ( J) =PYP1 ( J) -B* (PYP ( J) -CYP ( J) ) 

MY ( J) =PY1 ( J) -B* (PY ( J) -CY ( J) ) 

20 CONTINUE 
C 

C FUNCTION EVALUATIONS 
C 

DO 30, J=1,N 

DIFF=DELO+MX (J) -MX ( J~l) 

XMAS=P IDEN* ((PI/8.)* (DIFF**2 . ) +DI FF* (ALO+O . 5 *DELO+0 . 5* ( 1 . 
$ -PI)* (MX ( J) -MX ( J- 1) ) ) +0 . 5* (4 . -P I) * ( (0 . 5*DIFF 

$ +TE)**2.))+EM 

YMAS=XMAS 

TERM1=UCT+MYP ( J) -MYP (J-l) 

TERM2=UC+MXP (J) -MXP (J~l) 

C 

XPP ( J) = (DENSA* (SQRT ( (TERM1**2 . ) + (TERM2**2 . ) ) ) * (CL ( J) 

$ *TERM1+CD ( J) *TERM2) +FDX ( J) *SIGN (1 . ,MXP (J) ) 

$ +XK* (MX ( J+l) ~2 . *MX ( J) +MX ( J-l) ) ) / 

$ XMAS 

C 

YPP ( J) = (DENSA* (SQRT ( (TERM1**2 . ) + (TERM2**2 . ) ) ) * (CL ( J) 

$ *TERM2+CD(J)*TERM1)+FDY(J)*SIGN(1.,MYP(J)) 

$ +YK* (MY (J+l) -2. *MY (J) +MY (J-l) ) ) / 

$ YMAS 

30 CONTINUE 
C 

DO 40, J=1,N 
C 

C CORRECTOR 
C 

CXP 1 ( J) =C* (9 . *XP 1 ( J) -XP3 ( J) +D* (XPP ( J) +2 . *XPP 1 ( J) 

$ -XPP2(J))) 

CXI ( J) =C* (9 . *X(J, I) -X ( J , 1-2) +D* (CXP 1 ( J) +2 . *XP 1 ( J) 

$ -XP2(J))) 

CYPl ( J) =C* (9 . *YPl ( J) -YP3 ( J) +D* (YPP ( J) +2 . *YPP1 ( J) 

$ -YPP2 (J) ) ) 

CY1 ( J) =C* (9 . *Y ( J , I) -Y ( J , 1-2) +D* (CYP 1 ( J) +2 . *YP 1 ( J) 

$ -YP2(J))) 

C 

C FINAL VALUE - 73 - 
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XP ( J) *=CXPl ( J) +AB* (PXP1 ( J) -CXP1 (J) ) 

X ( J , 1+1) =CX1 ( J) +AB* (PX1 ( J) -CXI (J) ) 

YP ( J) =CYP1 ( J) +AB* (PYP1 ( J) -CYP1 ( J) ) 

Y ( J , 1+1) “CY1 ( J) +AB* (PYl (J) -CY1 (J) ) 

40 CONTINUE 
C 

DO 45, J*=0,N+1 
PX(J)“PX1(J) 

PXP(J)“PXP1(J) 

PY (J) “PYl (J) 

PYP (J) *=PYP1 (J) 

CX(J)“CX1(J) 

CXP(J) «CXP1(J) 

CY(J)-CY1(J) 

CYP (J) =CYP1 (J) 

XP4(J)=XP3(J) 

XP3(J) =XP2(J) 

XP2(J)=XP1(J) 

XP1(J)=XP (J) 

YP4(J)=YP3(J) 

YP3(J)=YP2(J) 

YP2(J)=YP1(J) 

YP1 (J)=YP (J) 

XPP4 ( J) =XPP3 (J) 

XPP3 ( J) =XPP2 (J) 

XPP2 (J) =XPP1 (J) 

XPP1(J) =XPP (J) 

YPP4 ( J) “YPP3 ( J) 

YPP3 ( J) =YPP2 (J) 

YPP2 ( J) “YPPl (J) 

YPP1 (J)=YPP (J) 

45 CONTINUE 
50 CONTINUE 
60 CONTINUE 
C 

RETURN 

END 

C 

C 

C 

SUBROUTINE STRAIN 
C 

C THIS SUBROUTINE COMPUTES THE STRAIN AMPLITUDE FOR 
C EACH OF THE CONVOLUTE TIPS. 

C 

COMMON/BLOK1/ PI , DELO , PIDEN , ALO , TE , EM, UC , UCT , N 
COMMON/BLOK3/ X,Y 
COMMON/BLOK5/ Al ,E,BOTX2,BOTY2 
COMMON/BLOK6/ SAMP 

DIMENSION X(0:6, 0:2048) ,Y(0:6, 0:2048) ,SAMP(5) , 
$ EXMAX (5) , EXMIN (5) , EYMAX (5) , EYMIN (5) , 

$ SMAX (5) , SMIN (5) 

C 

IF (N .GT. 1) GO TO 2 
AMULT=0. 1 
GO TO 4 

2 AMULT=0 . 13+0 . 29*FLOAT (N~2) 

4 BOTX2=AMULT*BOTX2 

BOTY2=AMULT*BOTY2 
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TERMX=E*A1*B0TX2 

TERMY=E*Al*B0TY2 

C 

DO 10 J=1 ,N 

EXMAX ( J) =X (J , 0) /TERMX 
EXMIN (J)=EXMAX(J) 

EYMAX ( J) =Y ( J , 0) /TERMY 
EYMIN ( J) “EYMAX ( J) 

10 CONTINUE 
C 

DO 30 I=1,20* Q 
DO 20 J “1 ,N 

EPX=X ( J , I) /TERMX 
EPY=Y ( J , I) /TERMY 
EXMAX ( J) =AMAX1 (EXMAX ( J) , EPX) 
EXMIN ( J) =AMIN1 (EXMIN ( J) , EPX) 
EYMAX ( J) =AMAX1 (EYMAX ( J) , EPY) 
EYMIN ( J) =AMINl (EYMIN ( J) , EPY) 
20 CONTINUE 
30 CONTINUE 
C 

DO 35 J=1,N 

SMAX (J) =EXMAX ( J) +EYMAX (J) 
SMIN ( J) =EXMIN ( J) +EYMIN ( J) 

35 CONTINUE 
C 

DO 40 J=1,N 

SAMP ( J) = (SMAX ( J) -SMIN ( J) ) /2 . 
40 CONTINUE 
C 

RETURN 

END 


SUBROUTINE FFT 

THIS SUBROUTINE TRANSFORMS THE SYSTEM RESPONSE 
PREDICTED BY THE MODEL INTO THE FREQUENCY DOMAIN 
SO THAT THE DOMINANT FREQUENCIES MAY BE DETERMINED. 

COMMON/BLOK1 / P I , DELO , P IDEN , ALO , TE , EM , UC , UCT , N 
COMMON / B LOK3 / X,Y 
COMMON/BLOK7/ XMAG , YMAG , OMEG 

DIMENSION X(0:6, 0:2048) ,Y(0:6, 0:2048) ,U(0:2048) , 
$ XMAG (0 : 2048) , YMAG (0 : 2048) , OMEG (0 : 2048) 

COMPLEX POWER , U , SUMX , SUMY , TERMX , TERMY , XT , YT 
C 

NUM=2048 
H=0. 00001 
T=H*NUM 
C 

DO 50 1=0, NUM 

OMEG (I) = (2 . *PI*FLOAT (i) ) /T 
50 CONTINUE 
C 

DO 1000 J=1,N 
C 

DO 100 1=0,1 

A=FLOAT (I) _ c 



POWER-CMPLX (0 . , (-2 . *P I*A) /FLOAT (NUM) ) 

U (I) =CEXP (POWER) 

100 CONTINUE 
C 

DO 200 1=2, NUM 
U(I)=U(I-1)*U(1) 

200 CONTINUE 
C 

DO 400 K=0,NUM-1 
SUMX=CMPLX (0 . ,0.) 

SUMY=CMPLX (0 . , 0 . ) 

L<0 300 M=0,NUM-1 

TERMX=X ( J , M) *U (MOD (M*K , NUM) ) 

TERMY=Y (J,M) *U (MOD (M*K, NUM) ) 

SUMX=SUMX+TERMX 
SUMY=SUMY+TERMY 
300 CONTINUE 

XT=H*SUMX 
YT=H*SUMY 
XMAG (K) =CABS (XT) 

YMAG (K) =CABS (YT) 

400 CONTINUE 
C 
C 

CALL FINDFRQ (J) 

C 

C 

1000 CONTINUE 
C 

RETURN 

END 

C 

C 

C 

SUBROUTINE FINDFRQ (J) 

C 

C THIS SUBROUTINE SEARCHES THROUGH THE FREQUENCY RESPONSE 
C DATA AND DETERMINES THE THREE MOST DOMINANT FREQUENCIES 
C IN BOTH THE LONGITUDINAL AND TRANSVERSE DIRECTIONS. 

C 

COMMON/BLOK1 / P I , DELO , P IDEN , ALO , TE , EM , UC , UCT , N 
COMMON /BLOK7 / XMAG , YMAG , OMEG 
COMMON / B LOK8 / XFRQ, YFRQ 

DIMENSION XFRQ (3, 5) , YFRQ (3, 5) , XMAG (0:2048) , YMAG (0 : 2048) , 
$ OMEG (0:2048) ,PK(3) 

C 

P=2.*PI 

C 

DO 10 K=l,3 
XFRQ (K, J) =0.0 
YFRQ (K, J) =0.0 
PK (K) =0.0 
10 CONTINUE 
C 

DO 30 1=2,60 

IF (XMAG(I) .LT. XMAG(I-l)) GO TO 30 
IF (XMAG (I) .LT. XMAG(I+1)) GO TO 30 
C 

IF (PK(1) .GT. XMAG(I)) GO TO 22 
PK (3)=PK(2) 
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PK(2)=PK(1) 

PK(1)=XMAG(I) 

XFRQ (3,J)=XFRQ(2,J) 

XFRQ (2 , J) =XFRQ (1 , J) 

XFRQ (1 , J) =OMEG (i) 

GO TO 30 

IF (PK (2) .GT. XMAG(I)) GO TO 24 
PK(3)“PK(2) 

PK(2) =XMAG(l) 

XFRQ (3, J)=XFRQ(2,J) 

XFRQ (2 , J) =OMEG (I) 

GO TO 30 

IF (PK (3) .GT. XMAG(I)) GO TO 30 
PK(3)=XMAG(I) 

XFRQ (3 , J) “OMEG (I) 

CONTINUE 

DO 35 K=1 ,3 
PK (K) =0 . 0 
CONTINUE 

DO 40 1=2,60 

IF (YMAG(I) .LT. YMAG(l-l)) GO TO 40 
IF (YMAG(I) .LT. YMAG(I+1)) GO TO 40 

IF (PK(1) .GT. YMAG(I)) GO TO 36 
PK(3)=PK(2) 

PK(2)=PK(1) 

PK(1)=YMAG(I) 

YFRQ (3 , J) =YFRQ (2 , J) 

YFRQ (2 , J) =YFRQ (1 , J) 

YFRQ (1 , J) =OMEG (I) 

GO TO 40 

IF (PK(2) .GT. YMAG(I)) GO TO 38 
PK(3) =PK (2) 

PK (2) =YMAG (i) 

YFRQ (3 , J) =YFRQ (2, J) 

YFRQ (2, J) =0MEG(I) 

GO TO 40 

IF (PK(3) .GT. YMAG(I)) GO TO 40 
PK(3)=YMAG(I) 

YFRQ (3, J)=0MEG(I) 

CONTINUE 

DO 45 K=1 ,3 

XFRQ (K, J) =XFRQ (K, J) /P 
YFRQ (K, J) =YFRQ (K, J) /P 
CONTINUE 

RETURN 

END 


SUBROUTINE CYCLE 



non 


C THIS SUBROUTINE CALCULATES THE NUMBER OF CYCLES 
C TO FAILURE FOR EACH CONVOLUTION. 

C 

COMMON/BLOK1/ PI , DELO , PIDEN , ALO , TE, EM, UC , UCT , N 
COMMON/BLOK5/ Al ,E,BOTX2,BOTY2 
COMMON/BLOK6/ SAMP 
COMMON/BLOK9/ CFAIL 
DIMENSION CFAIL (5) , SAMP (5) 

REAL NF 
C 

B— 0.1 
C— 0.7 

D-C12. *150000. )/E 
C 

DO 30 J=1 ,N 
NF=10. 

NUM=1 

10 FNF-D* ( (2 . *NF) **B) + ( (2 . *NF) **C) -SAMP ( J) /2 . 

FNF1-2 . *B*D* ( (2 . *NF) ** (B-l . ) ) +2 . *C* ( (2 . *NF) ** (C-l .) ) 
DELNF=FNF / FNF 1 

IF (ABS(DELNF) .LE. 0.5) GO TO 20 

NF=NF-DELNF 

NUM=NUM+1 

IF (NUM .GT. 50000) GO TO 40 
GO TO 10 

20 CFAIL (J)=NF 

30 CONTINUE 

GO TO 45 

40 PRINT*, ' CONVERGENCE PROBLEM IN CYCLE. ' 

45 RETURN 

END 


SUBROUTINE DAMAGE 
C 

C THIS SUBROUTINE COMPUTES THE DAMAGE DONE TO 
C EACH CONVOLUTION DURING SAMPLING PERIOD WHICH 
C INDICATES THE CONVOLUTION MOST LIKELY TO 
C FAIL FIRST. 

C 

COMMON/BLOK1 / P I , DELO , P IDEN , ALO , TE , EM , UC , UCT , N 

COMMON / B LOK8 / XFRQ.YFRQ 

COMMON/BLOK9/ CFAIL 

COMMON/BLOKIO/ NFAIL, FMAX, FMIN 

DIMENSION XFRQ (3,5), YFRQ (3,5) .CFAIL (5) ,D(5) 

C 

DO 20 J=1,N 
SUM=0 . 0 
DO 10 K=1 ,3 

TERM=XFRQ (K, J) +YFRQ (K, J) 

SUM=SUM+TERM 
10 CONTINUE 

SUM=0. 02048 *SUM 
D ( J) =SUM/ CFAIL ( J) 

20 CONTINUE 
C 

DTEST-0.0 

DO 30 J=1,N - 78 - 

IF (D(J) .LT. DTEST) GO TO 30 



non non 


dtest=d(j) 

NFAIL-J 
30 CONTINUE 
C 

XFMAX-0.0 
XFMIN= 1000000.0 
YFMAX=0 . 0 
YFMIN= 1000000.0 
DO 40 K=l,3 

IF (XFRQ (K, NFAIL) .EQ. 0.0) GO TO 40 
XFMAX=AMAX1 (XFRQ (K, NFAIL) ,XFMAX) 
XFMIN=AMIN1 (XFRQ (K, NFAIL) ,XFMIN) 

40 CONTINUE 

DO 50 K=l,3 

IF (YFRQ(K, NFAIL) .EQ. 0.0) GO TO 50 
YFMAX=AMAXl (YFRQ (K, NFAIL) , YFMAX) 
YFMIN=AMINl (YFRQ (K, NFAIL) , YFMIN) 

50 CONTINUE 
C 

FMAX=AMAXl (XFMAX, YFMAX) 

FMIN=AMINl (XFMIN, YFMIN) 

C 

RETURN 

END 


SUBROUTINE RESULT 

THIS SUBROUTINE PRODUCES THE DESIRED OUTPUT. 

COMMON / B LOK 1 / PI , DELO , PIDEN , ALO , TE , EM, UC , UCT , N 
COMMON/BLOK9/ CFAIL 
COMMON/BLOKIO/ NFAIL, FMAX.FMIN 
DIMENSION CFAIL (5) 

C 

TMAX=CFAIL (NFAIL) /FMIN 
TMIN=CFAIL (NFAIL) /FMAX 
C 


PRINT*, 

» ** 

NUMBER OF CYCLES 

TO FAILURE: 

PRINT*, 

f ** 

' .CFAIL (NFAIL) 


PRINT*, 

1 



PRINT*, 

t ** 

MINIMUM TIME TO 

FAILURE: ' 

PRINT*, 

» ** 

' , TMIN , ' SEC.' 


PRINT*, 

1 ** I 



PRINT*, 

1 ** 

MAXIMUM TIME TO 

FAILURE: ' 

PRINT*, 

I ** 

1 , TMAX , ' SEC.' 



C 

RETURN 

END 
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APPENDIX III 


SAMPLE DATA FROM PROGRAM 

The following charts show the cycles to failure, (Nf), taken from several 
test runs of the computer program. In each case the fluid flowing through the 
bellows is taken to be water with a density of 0.036 lbm/in 3 . N c is the 
number of convolutes and U c is the flow velocity in ft/s. 

Bellows 1 

Convolute pitch = 0.34 in. 

Convolute height = 0.492 in. 

Ply thickness = 0.015 in. 

Number of plys = 1 

Bellows inside diameter = 4.55 in. 

Elastic modulus = 29 x 10 8 psi 

Material density = 0.286 lbm/in 3 


Cycles to Failure 



2 

3 

4 

5 

10 

1.9 x 10 7 

1.3 x 10 9 

2.1 x 10 9 

5.9 x 10 9 

20 

1.6 x 10 6 

8.0 x 10 6 

4.5 x 10 5 

2.0 x 10 5 

30 

4.4 x 10 6 

1.4 x 10 8 

1.7 x 10 6 

4.1 x 10 5 

40 

4.7 x 10 7 

2.9 x 10 9 

3.5 x 10 7 

4.2 x 10 6 

50 

8.8 x 10 6 

6.0 x 10 7 

3.6 x 10 6 

8.6 x 10 6 

60 

2.4 x 10 5 

1.2 x 10 6 

6.2 x 10 5 

2.5 x 10 5 
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Bellows 2 


Convolute pitch = 0.3 in. 

Convolute height = 0.395 in. 

Ply thickness = 0.01 in. 

Number of plys = 2 
Bellows inside diameter = 4.47 in. 
Elastic modulus = 29 x 10® psi 
Material density = 0.286 lbm/in 8 


Cycles to Failure 


\«c 

u c 

2 

3 

4 

5 



10 

3.3x 10 8 

4.0 x 10 10 

8.9 x 10 10 

2.3 x 10 10 

20 

3.3 x 10 7 

3.4 x 10 8 

5.5 x 10 7 

8.3 x 10 6 

30 

9.7 x 10 7 

2.3 x 10 9 

6.1 x 10 7 

7.8 x 10 6 

40 

1.8 x 10 8 

5.0 x 10 10 

1.4 x 10 9 

2.2 x 10 8 

50 

2.8 x 10 7 

8.9 x 10 8 

1.8 x 10 8 

3.3 x 10 7 

60 

4.0 x 10 6 

1.6 x 10 7 

2.1 x 10 7 

8.4 x 10 6 











MSFC TEST BELLOWS 


BELLOWS 

NO. 

V F 

(FPS) 

TIME TO 
FAILURE 
(SEC) 

CYCLES TO 
FAILURE 
(X10 6 ) 

MATERIAL 

MSSm 

VARIATION 

5002-1 

62 

148-158 

.17-. 18 

321 


ELBOW 

5002-2 

72 

525-535 

.70-. 71 

321 


ELBOW 

5002-3 

48 

313-333 

.34-, 36 

321 

29,000 

ELBOW 

5002-4 

56 

2255-2309 

2. 5-2. 6 

321 

27,800 

ELBOW 

5002-5 

70 

10460-10585 

16.0-16.1 

321 

26,500 

ELBOW 

5002-6 

70 

730-736 

.95-. 96 

321 

28,100 

ELBOW 

5005-1 

33.2 

1841-1863 

.78-. 79 

321 

28 , 200 


5005-2 

65.2 

75-100 

— 

321 

*33 , 000 

ELBOW 

5006-10 

77.7 

1500-1530 

1.7 

321 

27,900 


5009-1 

40.5 

2728-2743 

1.8 

321 

27,900 


5011-1 

4 9.4 

3107-3133 

1.9 

321 

27,600 


5013-2 

50 

2007-2027 

1.3 

321 

28,000 


5013-3 

45 

749-760 

.47-. 48 

321 

28,800 


5028-1 

65.8 

787 


21-6-9 - 

*33,000 

ELBOW 

5028-2 

66.6 

150 

.16 

21-6-9 

37,000 

ELBOW 

5034-1 

82.5 

1672 


21-6-9 

31,200 

X 

5034-2 

84 

1427-1467 

■ItS " 

21-6-9 

31,200 


5034-3 

84.2 

3391-2402 

5.0 

21-6-9 

31,000 

% 

5034-4 

84.4 

3683 

5.3 

21-6-9 

31,000 

% 

5034-5 

65.5 

353-373 

. 37-.39 

21-6-9 

34,000 

-x 

5034-6 

85.5 

226-256 

.28-. 32 

21-6-9 

34,200 

X 

5034-7 

85 

2199 

2.9 

21-6-9 

31,000 

X 

5034-8 

85.8 

560-590 

.74-. 78 

21-6-9 

31,900 


5034-10 

70.9 

3580-3760 

3. 9-4.1 

21-6-9 

31,000 


5034-13 

79.5 

1375-1495 

1. 6-1.8 

21.6-9 

31,200 


5034-15 

69 

4180 

4.5 

21-6-9 

31,000 


AG01 

57.8 

1604-1634 

1.7 

21-6-9 

31,400 


AG#2 

74.8 

20099-20124 

«... 

21-6-9 

.... 


AG// 3 

80.5 

11248-11278 



21-6-9 

— 


5035-1 

113 

1976-1996 

4.4-4. 5 

21-6-9 

31,000 


5035-2 

110 

437-477 

.98-1.1 

21-6-9 

31,800 


5035-3 

94 

7916-8036 

16.3-16.6 

• 

21-6-9 

31,000 

X 


FROM REFERENCE 2. 
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